{
 "metadata": {
  "language": "Julia",
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "A simple example of fitting a nonlinear regression and nonlinear mixed-effects model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First we load the `NLreg` package and and the `DataFrames` package. The `readtable` function in creates a `DataFrame` from a comma-separated or tab-separated file's contents.\n",
      "The file can be compressed, as shown here."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "using DataFrames,NLreg\n",
      "const sd1 = readtable(Pkg.dir(\"NLreg\",\"data\",\"sd1.csv.gz\"));"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "Warning: redefining constant sd1\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Nonlinear regression models are represented as types in the `NLreg` package.  The type can implement several constructors.  Here we generate a `BolusSD1` model from a formula expression and a data frame."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m = BolusSD1(CONC ~ TIME, sd1);\n",
      "typeof(m)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "BolusSD1{Float64} (constructor with 1 method)"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This model is a concrete type inheriting from the abstract type `PLregModF`, a partially-linear regression model function."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "typeof(m) <: PLregModF"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "true"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- Add information here about the size of a model, once those methods are defined\n",
      "- Also add `formula` methods\n",
      "\n",
      "The formula indicates that the response is `CONC` and the first (and only) covariate is `TIME`. The names of the parameters for the model are provided by `pnames`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pnames(m)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "2-element Array{ASCIIString,1}:\n",
        " \"V\"\n",
        " \"K\""
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The first parameter, `V`, representing the volume of distribution, is conditionally linear. The second parameter, `K`, the rate constant is nonlinear."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Fitting a nonlinear regression model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A partially linear least squares fit is represented by the type `PLinearLS`, which is constructed from a nonlinear regression model and, optionally, initial values for the parameters."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If initial values for the parameters are not given the model's `initpars` method is used to generate initial values."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "show(initpars(m))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[0."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2793108696162248]"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that the initial estimates are for the nonlinear parameters only.\n",
      "\n",
      "The `fit` function fits the model using the Golub-Pereyra algorithm.  An optional second argument determines if verbose output is generated."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pl = fit(m,true)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Iteration:  1, rss = 111.039, cvg = 0.00367759 at [1.20139,0.279311]\n",
        "   Incr: [-0.0329797]  f = 1.0, rss = 110.597\n",
        "Iteration:  2, rss = 110.597, cvg = 1.53792e-6 at [1.14412,0.246331]\n",
        "   Incr: [-0.000604764]  f = 1.0, rss = 110.597\n",
        "Iteration:  3, rss = 110.597, cvg = 6.08026e-9 at [1.14303,0.245726]\n",
        "   Incr: [-3.79507e-5]  f = 1.0, rss = 110.597\n",
        "\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "Nonlinear least squares fit to 580 observations\n",
        "\n",
        "     Estimate Std.Error t value Pr(>|t|)\n",
        "V     1.14296 0.0495656 23.0595  < eps()\n",
        "K    0.245688 0.0202414 12.1379  < eps()\n",
        "\n",
        "Residual sum of squares at estimates: 110.597\n",
        "Residual standard error = 0.43743 on 578 degrees of freedom"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There are several extractor methods such as `coef`, `coeftable`, `deviance`, `pnames`, `stderr` and `vcov` for this type."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "coeftable(pl) # the coefficient table shown above"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "     Estimate Std.Error t value Pr(>|t|)\n",
        "V     1.14296 0.0495656 23.0595  < eps()\n",
        "K    0.245688 0.0202414 12.1379  < eps()\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "deviance(pl)  # residual sum-of-squares"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "110.5972863376179"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this model we used the elimination rate constant, `K`, as a parameter.  Suppose that instead we wish to use the logarithm of the rate constant as a parameter.  We could create a new model type or we could combine the existing model with a parameter transformation.\n",
      "\n",
      "- Consider how to code up parameter transformations of both parameters."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pl2 = fit([LogTr]*BolusSD1(CONC ~ TIME, sd1))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 19,
       "text": [
        "Nonlinear least squares fit to 580 observations\n",
        "\n",
        "        Estimate Std.Error  t value Pr(>|t|)\n",
        "V        1.14296 0.0495654  23.0595  < eps()\n",
        "log(K)   -1.4037 0.0823867 -17.0379  < eps()\n",
        "\n",
        "Residual sum of squares at estimates: 110.597\n",
        "Residual standard error = 0.43743 on 578 degrees of freedom"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The two fits are equivalent in that they produce the same sum of squared residuals and the parameter estimates can be transformed from one scale to the other."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "using Base.Test\n",
      "@test_approx_eq_eps coef(pl) [coef(pl2)[1], exp(coef(pl2)[2])] 1e-5"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Fitting a simple nonlinear mixed-effects model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A simple nonlinear mixed-effects model has a random effect for each parameter for each group (e.g. subject).  Uncorrelated random effects are indicated by specifying a `Diagonal` $\\Lambda$ matrix.  Correlated random effects user a `Triangular` $\\Lambda$ matrix.  The model form and initial estimates of the fixed-effects parameters can be taken from a fitted `NonlinearLS` object.\n",
      "\n",
      "- use an extractor function for mmf"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nm = fit(SimpleNLMM(pl2,array(sd1[:ID]),Diagonal))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "LoadError",
       "evalue": "type CompositePLModF has no field mmf\nwhile loading In[21], in expression starting on line 1",
       "output_type": "pyerr",
       "traceback": [
        "type CompositePLModF has no field mmf\nwhile loading In[21], in expression starting on line 1",
        " in updtmu! at /home/bates/.julia/NLreg/src/nonlinreg.jl:103",
        " in prss! at /home/bates/.julia/NLreg/src/simplenlmm.jl:54",
        " in pnls! at /home/bates/.julia/NLreg/src/simplenlmm.jl:100",
        " in pnls! at /home/bates/.julia/NLreg/src/simplenlmm.jl:117",
        " in fit at /home/bates/.julia/NLreg/src/nlmm.jl:51"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nm1 = fit(SimpleNLMM(pl2,array(sd1[:ID]),Triangular))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "LoadError",
       "evalue": "type Triangular is immutable\nwhile loading In[22], in expression starting on line 1",
       "output_type": "pyerr",
       "traceback": [
        "type Triangular is immutable\nwhile loading In[22], in expression starting on line 1",
        " in SimpleNLMM at /home/bates/.julia/NLreg/src/simplenlmm.jl:37"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Although not currently shown in the output, the within-subject correlation of the random effects for this model fit is `-0.091`.  A likelihood ratio test comparting `nm` versus `nm1` has a p-value of 39.6% from which we would conclude that the more complex model, `nm1`, does not provide a sufficiently better fit to justify the additional parameter."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}